// Copyright Mihai Preda

#include "Trig.h"

#include <cassert>

using namespace std;

namespace {

pair<u32, u32> splitTwos(u32 n) {
  int nTwos = 0;
  while (n % 2 == 0) {
    ++nTwos;
    n /= 2;
  }
  return {n, 1u << nTwos};
}


const pair<u32, u32> COS_TAB[8] = {{9, 0},  {3, 0},  {57, 3}, {21, 2}, {1, 0}, {13, 1}, {11, 1}, {19, 3}};

array<double, 8> COS[4] {
/* 9 */   {1,-0.0038077177473338575,2.4164524072268794e-06,-6.1341124777301022e-10,8.3417746250225797e-14,-7.0584717738912358e-18,4.0721382482643969e-22,-1.6919074352334798e-26,},
/* 143 */ {1,-1.5082651353809108e-05,3.7914395310092149e-11,-3.8123307049719434e-17,2.0535733802621082e-23,-6.8829557033353699e-30,1.5726316782026116e-36,-2.5487666383700158e-43,},
/* 147 */ {1,-1.4272994471472185e-05,3.395306186377832e-11,-3.2307457617727589e-17,1.6468720082881016e-23,-5.2235052014489976e-30,1.1294228181124433e-36,-1.7340690917912537e-43,},
/* 285 */ {1,-3.7971700527429051e-06,2.4030834015746363e-12,-6.0832775511763575e-19,8.2497283518972138e-26,-6.9612521485946193e-33,4.0052215406851532e-40,-1.6706346039780075e-47,},
};

array<double, 8> SIN[4] {
{0.087266462599716474,4.877106420722032e-18,-0.00011076201946266205,4.2175050723802824e-08,-7.6471756690855468e-12,8.0884113884997494e-16,-5.5994155932178749e-20,2.7063667699130723e-24,},
{0.005492294848933205,-2.0584448006122426e-20,-2.761278944626035e-08,4.1647407612360722e-14,-2.9912063260820776e-20,1.253203156624991e-26,-3.4364810453922909e-33,6.5792419764737163e-40,},
{0.0053428446489622331,2.4549457241226296e-19,-2.5419464045524189e-08,3.6281186978912718e-14,-2.4659103852004987e-20,9.776644053115473e-27,-2.5369946093312833e-33,4.5962859041262069e-40,},
{0.0027557830294647309,-2.4023522406397623e-20,-3.4880589304468587e-09,1.3244752912880034e-15,-2.3948847185853837e-22,2.5260507823548698e-29,-1.7438811621898675e-36,8.4054226417956712e-44,},
};

template<size_t N>
array<double, N> scaleSin(const array<double, N>& v, double f) {
  array<double, N> ret;
  ret.at(0) = v.at(0) * f;
  auto f2 = f * f;
  for (u32 i = 1; i < N; ++i) {
    ret.at(i) = v.at(i) * f;
    f *= f2;
  }
  return ret;
}

template<size_t N>
array<double, N> scaleCos(const array<double, N>& v, double f) {
  array<double, N> ret;
  auto f2 = f * f;
  f = 1;
  for (u32 i = 0; i < N; ++i) {
    ret.at(i) = v.at(i) * f;
    f *= f2;
  }
  return ret;
}

} // namespace

TrigCoefs trigCoefs(u32 n) {
  auto [mid, twos] = splitTwos(n);
  assert(twos >= 1 && (twos & (twos - 1)) == 0 && (twos % 4 == 0));
  assert(mid % 2 != 0 && 1 <= mid && mid <= 15);

  u32 pos = (mid - 1) / 2;

  assert(pos <= 7);
  auto [mul, idx] = COS_TAB[pos];

  return {mul, scaleSin(SIN[idx], 1.0/(twos / 4)), scaleCos(COS[idx], 1.0/(twos / 4))};
}
